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In  the  design  of  Fast  Multipole  Methods  (FMM)  for  the  numerical  solution  of 
scattering  problems,  a  crucial  step  is  the  diagonalization  of  translation  operators 
for  the  Helmholtz  equation.  These  operators  have  analytically  simple,  physically 
transparent,  and  numerically  stable  diagonal  forms.  It  has  been  observed  by  sev¬ 
eral  researchers  that  for  any  given  precision  e,  diagonal  forms  for  the  translation 
operators  for  the  Helmholtz  equation  are  not  unique,  and  that  some  choices  lead 
to  more  efficient  FMM  schemes  than  others.  As  is  well-known,  original  single- 
stage  FMM  algorithms  for  the  Helmholtz  equation  have  asymptotic  CPU  time 
requirements  of  order  0(n3/2),  where  n  is  the  number  of  nodes  in  the  discretiza¬ 
tion  of  the  boundary  of  the  scatterer;  two-stage  versions  have  CPU  time  estimates 
of  order  0(n4/3);  generally,  k- stage  versions  have  CPU  time  estimates  of  order 
0(n(fc+2U(fc+1)).  However,  there  exist  choices  of  diagonal  forms  leading  to  single- 
stage  FMM  algorithms  with  CPU  time  requirements  of  order  0(n4/3),  two-stage 
schemes  with  CPU  time  requirements  0(n5/4),  etc.  In  this  paper,  we  construct 
such  diagonal  forms  in  two  dimensions.  While  the  construction  of  this  paper  is  in 
no  sense  optimal,  it  is  rigorous  and  straightforward.  Our  numerical  experiments 
indicate  that  it  is  within  a  factor  of  two  of  being  optimal,  in  terms  of  the  number 
of  nodes  required  to  discretize  the  translation  operator  to  a  specified  precision  t. 
The  procedure  is  illustrated  with  several  numerical  examples. 
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Sparse  Diagonal  Forms  for  Translation  Operators  for  the  Helmholtz 

Equation  in  Two  Dimensions 


1  Introduction 

In  the  design  of  Fast  Multipole  Methods  (FMM)  for  the  numerical  solution  of  scattering  prob¬ 
lems,  a  crucial  step  is  the  diagonalization  of  translation  operators  for  the  Helmholtz  equation. 
These  operators  have  analytically  simple,  physically  transparent,  and  numerically  stable  diag¬ 
onal  forms.  Once  the  latter  are  constructed,  the  design  of  FMM  schemes  is  straightforward; 
the  simplest  “single-stage”  algorithms  have  CPU  time  requirements  of  order  0(n3/2),  where 
n  is  the  number  of  nodes  in  the  discretization  of  the  problem.  Two-stage  schemes  have  CPU 
time  requirements  of  order  0(n4/3);  generally,  fc— stage  schemes  have  CPU  time  requirements 
of  order  0(n(fc+2)/^+1)). 

It  has  been  observed  by  several  researchers  (see,  for  example,  [4,  2,  3])  that  for  any  given 
precision  e,  diagonal  forms  for  the  translation  operators  for  the  Helmholtz  equation  are  not 
unique,  and  that  some  choices  lead  to  more  efficient  FMM  schemes  than  others.  In  fact, 
there  exist  choices  of  diagonal  forms  leading  to  single-stage  FMM  algorithms  with  CPU  time 
requirements  of  order  0(n4/3),  two-stage  schemes  with  CPU  time  requirements  0(n5//4),  etc. 

Due  to  space  limitations,  we  will  not  describe  here  the  FMM  for  the  Helmholtz  equation, 
referring  the  reader  to  [6],  [7].  We  will  observe  that  the  functions  b^e  :  [0 ,  2tt ]  — ►  C  defined 
below  are  a  generalization  of  functions  vn  :  [0,2tt]  — ►  C  of  [6].  However,  while  the  functions  vn 
(see  (3.25)  in  [6])  are  nowhere  small  on  the  interval  [0, 2i r],  the  functions  are  negligibly  small 
on  most  of  their  interval  of  definition.  The  cost  of  this  sparsity  is  a  somewhat  higher  frequency 
content  of  the  functions  b£  £;  both  the  functions  u„,  b^£  are  trigonometric  polynomials  of  finite 
order,  and  the  order  of  6£e  is  1.5  times  higher  than  the  order  of  vn  under  similar  conditions. 
The  result  of  this  trade-off  is  a  reduction  in  the  cost  of  the  algorithm  (see  the  preceding  two 
paragraphs).  For  a  more  detailed  motivation  for  the  development  of  improved  translation 
operators  for  the  Helmholtz  equation,  we  refer  the  reader  to  the  papers  [4,  2,  3]. 

The  purpose  of  this  paper  is  to  construct  such  diagonal  forms  in  two  dimensions.  The  con¬ 
struction  is  intended  to  be  reasonably  rigorous;  it  is  also  quite  simple.  However,  the  proof  that 
the  resulting  translation  operator  is  sparse  is  quite  long  and  somewhat  technical.  Thus,  discus¬ 
sion  is  deliberately  conducted  on  two  levels.  First,  we  formulate  the  problem  (Subsection  1.1), 
and  describe  a  solution  (Subsection  1.2).  The  solution  of  Subsection  1.2  has  been  chosen  for 
analytical  simplicity,  rather  than  for  its  numerical  properties.  A  more  numerically  attractive 
solution  is  described  in  Subsection  5.1  and  illustrated  in  Subsection  5.1  by  numerical  examples. 
Thus,  a  reader  who  is  not  interested  in  the  proofs  may  want  to  read  Subsections  1.1,  1.2,  and 
turn  to  Subsections  5.2,  5.2,  possibly  after  reading  Subsection  1.3  (informal  description  of  the 
construction). 

Otherwise,  the  structure  of  the  paper  is  as  follows.  In  Section  2,  we  summarize  the  known 
facts  from  analysis  to  be  used  in  the  remainder  of  the  paper.  In  Section  3,  we  develop  the 
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requisite  analytical  apparatus.  In  Section  4,  we  prove  that  the  functions  introduced  in  this 
section  do  in  fact  satisfy  the  conditions  1.  -  3.  of  this  section.  Finally,  Section  5  contains  a 
slightly  modified  construction  of  the  functions  (quite  similar  to  that  of  this  section,  but 
leading  to  somewhat  faster  computations),  illustrated  by  several  numerical  examples. 

Remark  1.1  The  mathematical  techniques  used  in  this  paper  are  limited  to  elementary  anal¬ 
ysis;  however,  the  constructions  we  use  are  fairly  involved.  Thus,  the  proof  of  the  principal 
analytical  result  of  this  paper  (Theorem  1.1)  consists  of  a  fairly  long  sequence  of  definitions  and 
lemmas.  Most  of  the  latter  follow  immediately  from  the  preceding  ones  and  from  the  relevant 
definitions,  and  in  such  cases  the  proofs  are  omitted. 

1.1  Statement  of  the  problem 

In  agreement  with  standard  practice,  we  will  denote  by  Jm  the  Bessel  function  of  the  first  kind 
of  order  m,  by  Hm  the  Hankel  function  of  order  m,  and  by  Im  the  modified  Bessel  function  of 
order  m  (see,  for  example,  [1],  Chapter  9). 

Suppose  that  r,  />,  £  are  three  real  numbers,  such  that 

0  <  4  -r<p,  (1) 

and 

0  <  £  <  1/10.  (2) 

We  will  denote  by  v  the  smallest  of  all  positive  integer  numbers  such  that 

Jj(r)  <  £.  (3) 

for  all  j  >  u.  The  purpose  of  this  paper  is  to  construct  a  function  b^c  :  [0, 27r]  — *•  C,  satisfying 
the  following  conditions. 

1.  There  exist  a  positive  integer  A  (independent  of  p  and  r,  but  possibly  dependent  on  £ )  and 

complex  a_A.„+i,  a-A-H-2,- * *?  <*-i>  ao,  «i >  •••>  a a-i/-i>  such  that 

«,(«)=  £  «<•«'*'.  « 

JZZ-X-l/ 

for  all  6  e  [—7 r,;r]. 

2.  For  all  j  €  [—2  •  u,  2  •  */], 

\  aj  -  Hk(p)  \<  e.  (5) 

3.  There  exist  two  numbers  p,  q,  independent  of  r  and  p  (though  possibly  dependent  on  £), 
such  that 

I  KM  l<  £  <6) 
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(7) 


for  all  6  6  [ — tt,  7r]  such  that 

\0\  + 

r  p 

Remark  1.2  From  the  point  of  view  of  asymptotic  CPU  time  estimates,  it  is  sufficient  to 
construct  functions  bp£  satisfying  the  conditions  1.-3.  above.  In  terms  of  actual  computation 
times,  it  is  critical  that  the  coefficients  A,  p,  q  be  as  small  as  possible.  In  the  construction  of  the 
following  subsection,  A  =  3;  Theorem  1.1  formulated  in  the  following  subsection  (to  be  proved 
in  Section  4)  provides  values  p  =  8  •  \f2  •  log(l/e),  q  =  8.  The  actual  values  are  considerably 
smaller,  as  can  be  seen  from  the  numerical  examples  presented  in  Section  5. 

Remark  1.3  In  this  paper  we  construct  several  solutions  of  the  Helmholtz  equation  that  are 
negligibly  small  over  most  of  the  complex  plane,  without  being  equal  to  zero.  We  will  also 
be  dealing  with  restrictions  of  such  solutions  on  circles  and  lines  in  the  plane.  Abusing  the 
terminology  somewhat,  when  we  say  that  the  support  of  some  function  is  contained  by  some 
region,  we  mean  that  outside  that  region,  the  absolute  value  of  the  function  is  smaller  than  a 
preselected  £. 


1.2  Construction  of  the  functions  bpT  e 

In  this  subsection,  we  construct  functions  e  satisfying  the  conditions  1.  -  3.  of  Subsection  1.1. 
We  will  denote  by  p  the  smallest  integer  such  that 

p  >  (r-3  +  c  •  r~3)3,  (8) 


with 


(9) 


*  =  log(-£), 

and  by  u  the  real  number  defined  by  the  formula 


We  will  define  two  positive  integer  numbers  m,  n  via  the  formulae 


(10) 


(11) 


m  =  6  •  p 
5 

n  ~  2  '  P' 

respectively,  and  by  the  function  [ — 7r,  7r]  — >  C,  defined  by  the  formula 


f'A°) 


stn((n  +  |M)  .  eu-(co,(e)-i) 
si'n(|) 


(12) 

(13) 


(14) 
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(15) 


Finally,  for  each  p  >  4  •  r,  we  will  define  the  function  b(.e  :  [— 7r,  7r]  — »  C  by  the  formula 


m/2 

j=-m/2 

where  (/£,r)j  denotes  the  y  -  th  Fourier  coefficient  of  the  function  fCyT  ■  [— jt]  — ►  C. 

The  following  theorem  states  that  the  function  6^£  satisfies  the  conditions  1.  -  3.  of  the 
preceding  subsection,  and  provides  estimates  for  the  coefficients  p,  q  in  (7).  Its  proof  is  the 
principal  purpose  of  Sections  3,  4  of  this  paper. 

Theorem  1.1  Suppose  that  r,  p,  e  are  three  real  numbers  satisfying  the  inequalities  (1),  (2), 
and  the  function  6£e  is  defined  by  (15).  Then 

1.  For  all  j  €  [-2  •  u,  2  •  i/], 

I  (K-h  -  Bt(p)  l<  t.  (16) 


2.  I  K,.W  l<  £ 

for  all  0  6  [— 7r,7r]  such  that 


8 • V2  ■  log(\)  r 

e  > - +  8  • 

r  p 


(17) 


(18) 


1.3  Outline  of  the  proof  of  Theorem  1.1 

Put  informally,  Theorem  1.1  states  that  given  real  numbers  r,  p  such  that  4  •  r  <  p,  there  exists 
a  function  b  :  [— 7T,jr]  — *  €  such  that 

a.  b  is  a  trigonometric  polynomial  of  order  m,  with  m  ~  3  •  r. 

b.  The  first  4  •  r  coefficients  in  the  Fourier  series  of  b  are  defined  by  the  formula 

(b)j  =  Hj(p)  (19) 

for  all  j  such  that  |  j  |<  2  •  r. 

c.  b(0)  is  small  for  all  0  outside  a  small  neighborhood  of  the  point  0  =  0.  More  specifically,  the 
size  of  the  region  around  0  where  6(0)  >  s  may  depend  on  e,  but  has  to  be  of  the  order  r/p. 

In  this  formulation,  it  is  clear  that  Theorem  1.1  is  not  at  all  obvious,  except  when  p  ~  r2, 
or  greater.  Indeed,  in  this  case, 

jt»M  ~  \/(  £)  ■  S'-rv  (20) 

(see  (39)),  and  the  problem  of  finding  a  function  satisfying  the  conditions  a.  -  c.  of  this  section 
becomes  a  classical  problem  of  designing  a  low-pass  filter  that  is  (almost)  band-limited  in  both 
time  and  frequency  domains. 


4 


When  p  is  considerably  smaller  than  r2  (which  is  normally  the  case  in  situations  involving 
the  FMM),  such  simple  asymptotic  techniques  do  not  work.  In  this  regime,  Theorem  1.1  is 
a  consequence  of  detailed  analytical  properties  of  Hankel  functions,  and  its  proof  has  to  take 
these  into  account.  In  this  paper,  we  observe  that  (4)  can  be  rewritten  as 

m/2 

KP,0)=  E  7j  •  Hj{p)  •  e*'-7'9,  (21) 

j=-m/ 2 

with  the  condition  (5)  assuming  the  form 

I  7j  ~  1  l<  £•  (22) 

Now,  we  change  our  point  of  view,  interpret  the  pair  ( p,6 )  in  (21)  as  polar  coordinates  of  a 
point  in  R 2,  and  define  the  mapping  Q  :  R2  — ►  C  via  the  formula 

Q(x,y)  =  b(p,0),  (23) 

with  ( p,6 )  the  polar  coordinates  of  the  point  (x,y)  €  R2.  Clearly,  in  this  interpretation,  Q  is  a 
solution  of  the  Helmholtz  equation  (24)  satisfying  the  radiation  condition  (27).  Thus,  the  proof 
of  Theorem  1.1  has  been  reduced  to  constructing  a  solution  to  the  equation  (24)  possessing 
certain  properties.  Once  such  a  solution  is  found,  the  function  b^e(0)  is  obtained  as  a  restriction 
of  Q  on  the  circle  of  radius  p. 

The  conditions  to  be  satisfied  by  Q  follow  immediately  from  the  conditions  a.  -  c.  above.  In 
addition  to  (22),  Q  must  look  like  a  beam  to  satisfy  (17).  Fortunately,  beam-like  solutions  of  the 
Helmholtz  equation  are  well-known;  a  typical  example  are  the  so-called  Gaussian  beams  (see 
Section  3  below).  In  this  paper,  we  obtain  functions  (21)  as  linear  combinations  of  Gaussian 
beams  (see  Section  4,  where  the  resulting  function  R2  — ►  C  is  denoted  by  Q™). 

The  last  problem  we  encounter  is  the  fact  that  Gaussian  beams  are  not  sufficiently  sharp  to 
satisfy  the  condition  c.  of  Subsection  1.1;  put  differently,  a  Gaussian  beam  that  is  sufficiently 
sharp  to  satisfy  the  condition  c.  is  singular  on  a  region  too  large  for  the  condition  a.  to  be 
satisfied.  Fortunately,  Gaussian  beams  can  be  modified  to  reduce  the  size  of  the  singular  region 
dramatically,  leaving  the  beam  almost  intact  away  from  the  singularity.  Section  3  is  largely 
devoted  to  this  construction,  which  is  referred  to  as  Modified  Gaussian  Beam. 

2  Analytical  Preliminaries 

In  this  section,  we  summarize  several  facts  from  analysis  to  be  used  in  the  sections  below.  All 
of  these  facts  are  either  well-known,  or  follow  immediately  from  well-known  facts. 

2.1  Notation. 


For  the  Helmholtz  equation 

VV  +  k2<f>  =  0 
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(24) 


we  will  define  the  potential  <j>^0  :  R2  \  {xo}  -*  C  of  a  unit  charge  located  at  the  point  xo  £  R2 
by  the  formula 

<(x)  =  #o(*||z-xo||),  (25) 

where  Ho  denotes  the  Hankel  function  of  order  zero.  We  will  define  the  potential  <t>XOih  a 
unity  dipole  located  at  x0  and  oriented  in  the  direction  h  £  R2  by  the  formula 

#.*(*)  -  -g,  Wl»  -  *.11)  •  (26) 

where  H\  denotes  the  Hankel  function  of  order  one.  In  most  cases,  a  potential  <f>  satisfying  the 
equation  (24)  in  an  unbounded  region,  also  satisfies  the  radiation  condition  at  oo,  i.e.  for  any 
x  €  R2,  there  exists  c  €  C  such  that 

Um  ■  x)  ■  e~iktW  ■  Vi  =  c,  (27) 

and  will  refer  to  functions  satisfying  the  equation  (24)  (and  in  unbounded  regions  -  also  the 
condition  (27)  )  as  radiation  potentials. 

Remark  2.1  In  the  remainder  of  the  paper,  we  will  be  assuming  that  the  Helmholtz  coefficient 
k  in  (24)  -  (27)  equal  to  1,  unless  explicitly  stated  otherwise. 

For  an  arbitrary  set  D  £  R2  and  a  point  x  £  R2,  we  will  denote  by  TX(D)  the  set  of  all 
points  y  in  R2  such  that  y  -  x  £  D. 

Given  a  set  D  C  R2  and  a  positive  real  number  r,  we  will  denote  by  Sr(D )  the  set  of  all 
points  x  £  R2  such  that  x  =  x0  +  y,  with  x0  €  D ,  and  y  some  vector  in  R2  such  that  ||y||  <  r. 
The  following  obvious  lemma  provides  a  bound  on  the  radius  of  ST{D)  given  the  radius  of  D. 

Lemma  2.1  Suppose  that  D  is  a  subset  of  R2,  and p  >  0  is  a  real  number,  such  that  ||xo||  <  p 
for  all  xq  £  D.  Then  ||x||  <  r  +  p  for  any  x  £  Sr(D ) 


In  Subsection  2.2,  we  will  need  the  following  lemma;  its  proof  is  an  exercise  in  elementary 
calculus,  and  is  omitted. 


Lemma  2.2  For  any  positive  real  x,t  and  natural  k, 


U  +  £)* 


<  e  *  e  2  *  . 


(28) 
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2.2  Elementary  properties  of  Bessel  functions 

As  is  well-known,  there  exist  two  functions  a,  (3  :  C  — ►  C,  such  that 

H0{z)  =  a(z)  +  0(z)  •  log(z)  (29) 

for  all  z  €  C;  the  functions  Ho,  Hi  are  connected  by  the  formula 

JflW  =  -£■ Betz).  (30) 

The  following  lemma  provides  a  crude  (but  sufficient  for  our  purposes)  estimate  of  the  absolute 
values  of  the  Hankel  functions  Ho ,  H\.  It  is  an  immediate  consequence  of  9.1.12,  9.1.13,  9.2.1, 
9.2.7  in  [1]. 

Lemma  2.3  For  any  z  6  C, 

p-Im(z) 

I  Ho(z)  |<  (31) 

VI*  I 

and 

+  (*» 

As  is  well  known  (see,  for  example,  [8]),  Jm  are  analytic  on  the  whole  complex  plane  for 
all  integer  values  of  m,  while  Hm  have  a  branch  cut  along  the  negative  real  axis,  and  become 
infinite  at  the  origin.  The  asymptotic  behavior  of  the  functions  Jm ,  Hm  for  large  m  is  given 
by  the  formulae 


hm  Jm(z)-(^)m- v/(2’rm)=  1, 

m— ►co  €2  v 

(33) 

lim  *„(*)•  (f 

m^oo  mv  '  v2 mJ  y/2 

(34) 

(see  [1],  9.3.1,  9.3.2,  9.1.3).  It  is  immediately  clear  from  (33)  that  the  functions  Jm(z)  decay 
rapidly  when  z  is  fixed  and  m  is  large.  However,  (33)  is  an  asymptotic  statement,  understating 
the  actual  rate  of  decay  of  (33)  when  m  is  only  slightly  greater  than  |  z  |.  For  purely  imaginary 
z,  a  dramatically  stronger  estimate  is  given  by  Lemma  2.7  below;  for  purely  real  z,  a  fairly 
tight  estimate  is  provided  by  the  following  lemma,  which  can  be  found  in  [8],  pp.  227,  255. 


Lemma  2.4  For  any  real  0  <  x  <  1  and  v  >  0, 


Jviy  •  a;)  < 


xv  •  e"'V/(1-12) 

(2  •  7T  •  •  (1  -  X 2)7  •  (1  +  y/(l  -  X2)Y 


(35) 


The  following  lemma  provides  a  simplified  version  of  (35).  Given  (35),  its  proof  is  an  exercise 
in  elementary  calculus,  and  is  omitted. 


7 


Lemma  2.5 

For  any  real  r  >  10  and  0  <  e  <  0.1, 

Jn(r)  <  e 

(36) 

for  any 

,1  1.0 
n  >  (n  +  c  •  r  3)°, 

(37) 

with  c,6  defined  by  (9),  (10). 


Remark  2.2  Obviously,  if  n  satisfies  the  inequality  (37)  and  v  is  defined  by  (3),  then  v  <n. 
For  large  z  and  fixed  m,  the  asymptotic  behavior  of  Jm(z),  Hm(z)  is  given  by  the  formulae 

=  (38) 

=  (39) 

when  z  — ►  oo,  as  long  as  Im(z )  >  0  (see  [1],  9.2.5,  9.2.7). 

We  will  need  the  behavior  of  Bessel  functions  in  one  more  asymptotic  regime,  as  provided 
by  the  following  lemma,  which  can  be  found  (in  a  slightly  different  form)  in  [8]. 


Lemma  2.6  For  any  integer  n  >  0, 

1  ra>  . .  nil  , 

2  2*/3 . 31/3  .  x  .  ni/3  *'  1  22/3  ■  31/3  ■  ir  •  n'/3  ' 

Furthermore , 

22/3  •  31/3  •  7r  •  n1/3  T,,  1 

hm  - — tt - J„(n)  =  1. 

n— *°°  r(|) 


(40) 


(41) 


As  is  well-known,  the  modified  Bessel  functions  Im  are  defined  by  the  formula 

Im(z)  =  i~n  ■  Jm(i  ■  Z)  (42) 

for  all  complex  2;  we  will  need  the  classical  formula 

OO 

ef(‘+D=  £  tk-Ik(z),  (43) 

k=  —  OG 

valid  for  all  pairs  2,  t  such  that  t  ^  0. 

It  is  well-known  that  once  n  >  2,  the  functions  Jn(z)  decay  rapidly  with  n  (see  (33),  (35)), 
for  all  complex  2.  What  appears  to  be  less  well-known,  is  that  when  2  is  purely  imaginary,  the 
decay  starts  at  n  ~  y/(2  ■  2).  The  following  lemma  provides  a  somewhat  crude  description  of 
the  behavior  of  In{x)  in  the  regime  y/{2  •  x)  <  n  <  x.  We  present  an  outline  of  the  proof  for 
this  lemma,  since  the  author  has  failed  to  find  it  in  the  literature. 
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(44) 


The  following  technical  lemma  is  obtained  from  the  preceding  one  by  elementary  algebraic 
manipulation. 

Lemma  2.8  Suppose  that  r,u,6  are  three  positive  real  numbers ,  and  n  is  an  integer  num¬ 
ber .  Suppose  further  that  6,  c,u  are  defined  by  (10),  (9),  (11),  respectively,  and  n  satisfies  the 
inequality  (37).  Then 

In(u)<eu-e~s.  (47) 

Finally,  we  will  need  two  well-known  integral  expressions  for  Bessel  functions,  given  by  the 
following  lemma. 

Lemma  2.9  For  any  integer  n  and  complex  z, 

j—n  r2n  ,  . 

Jn(Z)  =i~-  €i-z-cos{6)  ,  »  (48) 

2  •  7T  Jq 

and 

In(z)  =  -L  .  [2*  ez-cos{e)  •  eitn  ede.  (49) 

2  •  7T  Jo 
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2.3  Green’s  formula  for  the  Helmholtz  equation 

The  following  theorem  is  a  special  case  of  the  famous  Green’s  formula.  It  can  be  found  (for 
example)  in  [5]. 

Theorem  2.10  Suppose  that  the  function  (p  :  R2  — ►  C  satisfies  the  Helmholtz  equation  (24) 
outside  the  region  Q,  with  boundary  F.  Suppose  further  that  it  satisfies  the  outgoing  radiation 
condition  (26)  at  oo.  Than  for  any  x  G  R2  \  ft, 

«*)  =  -j  ■  fm  ■  If  (<,*)  +  (50) 

with 

G(x,y)  =  H0(k-  \\x-  y||)  (51) 

for  any  x,y  €  R2  such  that  x  ^  y;  the  integration  in  (50)  is  with  respect  to  the  arc  length. 

2.4  Partial  wave  expansions  of  radiation  potentials 

Suppose  that  a  function  ip  :  R2  — >  C  satisfies  the  Helmholtz  equation  (24)  outside  the  disk 
D  of  radius  R  with  the  center  at  the  point  xq  €  R2,  and  that  it  also  satisfies  the  radiation 
condition  (26)  at  oo.  Then  there  exists  a  unique  sequence  a  =  {am},  m  —  0, 1, 2,  *  •  *,  such  that 
for  any  x  6  R2  \  D, 

+oo 

tf(*)=  £  c*m-Hm(kp).eime.  (52) 

m  —  —  co 

In  the  above  formula,  p  =  ||x  —  xo||  and  6  is  the  angle  between  the  vector  x  —  xq  and  the  x  axis. 

A  derivation  of  the  formula  (52)  can  be  found,  for  example,  in  [5j;  we  will  refer  to  expansions 
of  the  form  (52)  as  if -expansions,  and  to  the  point  zo  as  the  center  of  the  expansion  (52). 

The  following  lemma  is  a  direct  consequence  of  the  formulae  (33),  (34).  It  establishes  the 
convergence  rate  of  the  expansion  (52). 

Lemma  2.11  If  D2  C  D  is  a  disk  of  radius  R2  >  R  with  the  center  at  xo  then  there  exists 
c  >  0  such  that  for  any  x  €  R2  \  &2  and  A  >  \k\  •  R, 

W*)-  £  Rm-Hm(kp)-eime\<c-(^-f.  (53) 

m=-N  2 

Remark  2.3  In  numerical  calculations,  the  expansion  (52)  is  truncated  after  a  finite  number 
of  terms,  and  the  resulting  expression  is  viewed  as  an  approximation  to  the  potential  ip.  As  is 
well-known,  if  we  want  to  approximate  ip  by  an  expansion  of  the  form  (53)  with  accuracy  £, 
we  have  to  choose 

N~R'\k\,  (54) 

i.e.  the  number  of  terms  in  the  approximation  is  almost  independent  of  e,  and  must  be  roughly 
equal  to  |Ar|  •  R. 
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2.5  Far-field  representations  of  radiation  potentials 

In  this  subsection,  we  introduce  an  alternative  form  of  the  expansion  (52),  possessing  a  simple 
physical  interpretation  simplifying  many  calculations  with  radiation  potentials. 

For  the  expansion  (52),  we  will  consider  a  function  FXo(ip )  :  [ — 7r,  7r]  —*•  C1  defined  by  the 
formula 

-F*o  WW  =  lim  ip(t  •  x  +  x0)  •  y/t  •  e~xkt  •  •  e* '*  (55) 

with  x  =  (cos  0,  sin  0).  Substituting  (39),  (52)  into  (55),  we  immediately  obtain 

+oo 

FX0(i>m=  E  Pme-^e™9,  (56) 

m~—oc 

which  provides  an  explicit  expression  for  FXo(ip)  v^a  the  coefficients  {/?m}.  Clearly,  (56)  defines 
a  unitary  mapping  connecting  the  coefficients  {c^}  in  the  expansion  (52)  with  the  with  the 
function  FX0  (0),  and  we  will  refer  to  FXo(ip)  as  the  far-field  representation  of  the  radiation 
potential  ip  with  the  origin  at  Xq . 

Obviously,  given  a  radiation  potential  (52),  its  far-field  representation  (55)  depends  on  the 
origin  xo .  The  following  lemma  describes  the  dependence;  its  proof  can  be  found  (for  example) 
in  [6]. 

Lemma  2.12  Suppose  that  the  radiation  potential  ip  is  defined  by  the  formula  (52),  and  xo, 
Xi  are  two  arbitrary  points  in  R2.  Suppose  further  that  FXo(ip),  FXl(ip)  are  the  far-field  repre¬ 
sentations  of  ip  with  origins  xq,  xj  respectively.  Then  for  any  6  €  [— 7r,7r], 

Fxime)  =  FXo{iP){0)  ■  (57) 

with  the  vector  w  £  R2  defined  by  the  formula 

w  =  (cos(O),  sin(6)).  (58) 

Remark  2.4  In  fact,  both  the  existence  of  asymptotic  representations  of  radiation  fields  and 
the  above  lemma  are  an  immediate  consequence  of  the  radiation  condition  (27).  A  detailed 
investigation  of  such  issues  can  be  found  in  [6]  in  the  two-dimensional  case,  and  in  [7]  in  the 
three-dimensional  one. 

The  following  lemma  is  an  immediate  consequence  of  the  formulae  (56),  (52).  It  provides 
an  explicit  formula  for  the  evaluation  of  a  radiation  potential  at  a  point,  given  its  far-field 
representation. 

Lemma  2.13  Suppose  that  a  function  ip  :  R2  — ►  C  satisfies  the  Helmholtz  equation  (24)  outside 
the  disk  D  of  radius  R  with  the  center  at  the  origin,  and  that  ip  also  satisfies  the  radiation 


11 


condition  (26)  at  oo.  Suppose  further  that  o  :  [ — 7r,  7r]  — >  C  is  the  far-field  representation  of  ip. 
Then  for  any  x  €  R2  \  D, 

+00 

rp(x)=  Y,  (59) 

m=-oo 

where  (p,  6)  are  the  polar  coordinates  of  x,  and  (5)m  denotes  the  m-th  term  of  the  Fourier  series 
of  o. 

2.6  Gaussian  Beams 

Gaussian  beams  are  solutions  of  the  Helmholtz  equation  that  are  obtained  as  potentials  of 
charges  with  complex  coordinates.  Such  a  potential  is  small  everywhere  outside  a  region  in  the 
plane  that  looks  like  a  spreading  beam,  and  inside  that  region  the  graph  of  the  absolute  value  of 
such  a  function  looks  like  the  normal  distribution;  hence  the  term  “Gaussian  beam”.  Because 
of  their  localized  nature,  Gaussian  beams  are  used  as  building  blocks  for  the  construction  of 
other  solutions  for  the  Helmholtz  equation. 

Suppose  that  u  >  0  is  a  real  number.  We  will  define  the  function  Gu  •  R2  — ♦  C  by  the 
formula 

Gu(x,y)  =  H0(/((x-i-u)2  +  y2)))-e-'t,  (60) 

and  refer  to  Gu  as  a  Gaussian  beam  with  base  ti,  oriented  horizontally  and  pointing  to  the 
right.  For  any  x0  6  R2,  we  will  denote  by  G„°  the  mapping  R2  -*  C,  defined  by  the  formula 

Gxu°(x)  =  Gu(x  +  x  0).  (61) 

Obviously,  G°  =  Gu,  and  Gx°  will  also  be  referred  to  as  a  Gaussian  Beam. 

The  following  lemma  provides  an  explicit  expression  for  the  far-field  representation  for  the 
Gaussian  beam  (60).  While  its  proof  is  very  simple,  we  provide  it,  since  it  reveals  important 
features  of  the  behavior  of  Gaussian  beams  at  large  distances  from  the  origin. 

Lemma  2.14  Suppose  that  u  is  a  positive  real  number.  Then  the  Gaussian  beam  (60)  is  a 
radiation  potential  outside  the  subset  Bu  of  R2,  consisting  of  all  pairs  (x,y)  such  that 

x  =0, 

—u  <  y  <  u. 

Furthermore,  the  far-field  representation  of  (60)  is  given  by  the  formula 

Fq{Gu)(8)  =  (64) 

for  all  6  €  [-ir,  it]. 


(62) 

(63) 
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Proof.  We  start  with  observing  that  the  equation 
(x  —  i  •  u)2  -f  y2  =  0 


(65) 


has  exactly  two  solutions: 


(  x  =  0,  y  =  -u),  (66) 

(  x  =0,  y  =  u).  (67) 


Thus,  the  function  Gu  :  R2  — ♦  C  has  logarithmic  singularities  at  the  points  (66),  (67),  connected 
by  a  branch  cut  (see  (29));  obviously,  outside  this  branch  cut  the  function  (60)  is  analytic  and 
satisfies  the  equation  (24). 

Suppose  now  that  x2  +  y2  >>  u2.  Using  the  Taylor  theorem,  we  have 


/((*  - i  • u )2  +  v2)  = 

/(x2  +  y2)./(l  2-i-X'U 

J(x2  +  y2)  -  t 


U* 


x 2  +  y2  x2  +  y 2 
x  •  u  1 


)  = 


+  0(- 


V(x2  +  y2)  W(*2  +  y2) 


:)  = 


r  -  i  •  u  -  cos(0 )  +  O(-), 
r 


(68) 


with  (r,  6)  the  polar  coordinates  of  the  point  (x,  y)  €  R2.  Now,  (64)  follows  from  the  combina¬ 
tion  of  (68),  (60),  (55),  (39). 

□ 


Remark  2.5  An  exercise  in  elementary  calculus  shows  that  for  any  u  >  0, 

e*.(eM(*)-i)  <  _e-*4  (69) 

for  all  6  €  [ — 7r,  ir],  and 

max  I  e«-(co.(tf)-i)  _  e-u4  .<  I  (70) 

eel-1 r,7r]  1  u  v  ; 

Both  bounds  (70),  (69)  are  quite  crude,  but  sufficient  for  our  purposes. 

3  Detailed  Analysis  of  Gaussian  Beams 

In  this  section,  we  develop  the  analytical  apparatus  to  be  used  in  Section  4  to  prove  Theorem  1.1. 
The  principal  tool  we  use  consists  of  the  well-known  Gaussian  beams;  however,  the  analysis  we 
use  is  somewhat  more  detailed  than  what  appears  to  be  present  in  the  literature. 
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3.1  Four  elementary  lemmas 

In  this  and  the  following  subsections,  we  analyze  the  spatial  structure  of  functions  of  the 
form  (60)  in  some  detail,  in  the  process  justifying  the  use  of  the  term  “Gaussian  beams”.  Lem¬ 
mas  3.1  -  3.4  below  provide  the  necessary  analysis.  Their  proofs  are  an  exercise  in  elementary 
calculus,  and  are  omitted.  We  start  with  several  additional  definitions. 

For  any  positive  real  u,  we  will  define  the  function  /  :  R2  — ►  C  by  the  formula 

/(*>  y)  =  \f((x  ~  i  •  «)2  +  y2)-  (71) 

For  any  positive  real  u,(3  such  that  (3  <  u,  we  will  define  two  regions  A,  B  in  R2,  as  follows. 

1.  The  region  AUt/s  consists  of  all  pairs  ( x,y ),  such  that 


y  l> 


V((x  2  +  /?2)-(u2-/J2)) 


2.  The  region  Bu>p  consists  of  all  pairs  (x,  y),  such  that 

^/T(»a  +  is2)  •  («3  -  /?2)) 


(see  Figure  1). 


(72) 


(73) 


Lemma  3.1  For  any  positive  u,/3  such  that  0  <  u, 

R2  =  Au>p\jBUt0.  (74) 

Furthermore,  for  any  (x,  y)  €  Au$, 

Im(f(x,y))<  /?,  (75) 

and  for  any  (x,  y )  €  BU}p, 

Jm(f(x,y))  >  (3,  (76) 

The  preceding  lemma  describes  precisely  the  part  AUtp,  of  the  plane  where  the  inequal¬ 
ity  (75)  is  satisfied;  the  following  one  provides  a  simplified  approximate  description  of  the 
region  AUip. 

Lemma  3.2  Suppose  that  under  the  conditions  of  the  preceding  lemma,  6  is  a  positive  real 
number,  such  that 

(3  =  u  -  6,  (77) 

and 

«  +  l<|.  (78) 
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Then 
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Lemma  3.4  Suppose  that  the  positive  real  numbers  x,  u,  6  are  such 

c  U 

6<r 


and 


x  >  2  •  u. 

Then 

Im(fu(x,  y))  <  u  -  6 
for  any  y  such  that 

I  y  I  >  4-V? 

X  y/u 

Proof.  Introducing  the  notation 
x 

-  -  V, 
u 

we  observe  that,  due  to  (88),  n  >  2,  and,  therefore, 

S{1  tOp.  <  V2. 

Now,  substituting  (91)  into  (90),  we  have 

|  y  |  4  •  y/S  2  ♦  y/2  •  y/6  •  u  » y/{l  -f  p1)  _ 

x  ^  ^  p  •  u3/2 

2  •  y/2  •  \/£  •  u  •  +  A*2  1  ^2)  _ 

p  *  u3/2 

2  *  \/2  •  •  u  •  +  x2) 

p  •  u3/2 

Multiplying  both  sides  of  (93)  by  z  and  using  (91)  again,  we  get 

2  •  \/2  •  \/J  ■  u  •  +  ^2  *  u2) 

l!,|> - -pjn - I  = 

2  •  \/2  •  \fl  •  \[{u 2  +  z2) 

Now,  (89)  follows  from  the  combination  of  (94)  and  Lemma  3.2. 


3.2  Geometry  of  Gaussian  beams 

In  this  subsection,  we  use  the  elementary  lemmas  of  the  preceding  one  to  describe  in  some 
detail  the  spatial  behavior  of  Gaussian  beams.  The  theorem  below  follows  immediately  from 
the  combination  of  Lemmas  3.3,  3.4  above  and  Lemma  2.3. 

Theorem  3.5  Suppose  that  the  positive  real  numbers  u,  6  are  such  that 

*<!•  (95) 

Then 

\Gu(x,y)\<e~s  (96) 

in  either  of  the  following  two  (intersecting)  regions: 

1.  x  €  [2  •  y/u,  2  •  u],  |  y  |>  2  •  i/(lO  •  S)  ■  \fu. 

2.  x  e  [2 -11,00],  \  y\>A$'X. 

Observation  3.1  The  above  theorem  has  a  very  transparent  physical  interpretation .  Specifi¬ 
cally,  a  Gaussian  beam  (60)  begins  to  look  like  a  beam  once  x  >  y/u.  While  x  £  [y/u,u],  the 
Gaussian  beam  virtually  does  not  expand .  At  approximately  x  =  uf  the  beam  begins  to  expand 
with  the  angle  of  expansion  roughly  4  •  y/S/yfu,  with  e~6  relative  error  of  our  measurements  (or 
calculations).  This  behavior  is  quite  obvious  in  Figure  3. 

4  Modified  and  Modulated  Gaussian  Beams 

4.1  Modified  Gaussian  beams 

According  to  Lemma  2.14,  the  Gaussian  beam  (60)  has  a  localized  far-field  representation. 
Specifically,  for  any  £  >  0, 

I  Fo(Gu)(0)  |<  £  (97) 

for  all  6  such  that 

|9|>V^|-l°g(l))  (98) 

(see  (69)).  On  the  other  hand,  Gu  has  logarithmic  singularities  at  the  points  (66),  (67),  and 
a  branch  cut  connecting  them.  In  other  words,  to  a  specified  precision,  the  support  of  the 
far-field  representation  of  a  Gaussian  beam  is  proportional  to  1  /\fu,  with  u  the  size  of  the 
region  where  the  beam  (60)  is  discontinuous.  This  very  large  region  of  discontinuity  turns  out 
to  be  a  major  problem  when  Gaussian  beams  are  used  as  bricks  for  the  construction  of  other 
solutions  of  the  equation  (24).  Fortunately,  there  exist  solutions  of  the  equation  (24)  almost 
exactly  coinciding  with  (60)  away  from  the  branch  cut  (62),  (63),  and  singular  on  a  region  of 
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size  roughly  *Ju.  We  will  refer  to  such  solutions  as  Modified  Gaussian  Beams  (MGBs),  and 
use  use  them  as  building  blocks  for  the  construction  of  solutions  of  the  equation  (24).  In  this 
subsection,  we  construct  the  MGBs,  and  prove  some  of  their  properties.  We  start  with  several 
definitions. 

For  an  arbitrary  pair  of  positive  real  numbers  u,  S,  we  will  denote  by  R\  the  rectangle  in 


the  plane  defined  by  the  four  vertices 

(l,-(ti+  l)),(l,(ti+  l)),(-l,(tt+  l)),(-l,-(u+  1)),  (99) 

by  R4  the  rectangle  defined  by  its  vertices 

(2,-(u  +  2)),(2,(u+  2)),(-2,  — (u  +  2)),(-2,  (u+  2)),  (100) 

by  the  rectangle  defined  by  the  vertices 

(1,-7),  (1,7),  (-1,7),  (-1,-7),  (10!) 

with 

7  =  2*  \J{2  •«•(£+  1)),  (102) 

and  by  R3  the  difference  R\  \  R2,  observing  that  R3  consists  of  two  rectangles,  with  vertices 
(1,-7), (!,-(«+  l)),(-l,-(«+  !)),(-!, -7),  (103) 

(1,7),  (!,(«+  1)), (-!,(«  +  !)),(-!, 7),  (104) 


respectively  (see  Figure  2).  We  will  denote  by  Ti,  T2,  T3,  T4  the  boundaries  of  the  regions  R\, 
f?2>  -R3,  Rii  respectively.  Whenever  a  function  is  to  be  integrated  over  one  of  these  boundaries, 
the  integration  will  always  be  assumed  to  be  with  respect  to  the  arc  length. 

We  will  denote  by  GUis  the  function  R2  \  R2  -*  C  defined  by  the  formula 

■  J  (G,(«)  •  ||(I,x)+  •  G(t,x))il,  (105) 

and  refer  to  GUys  as  a  Modified  Gaussian  Beam.  For  any  x0  £  R2,  we  will  denote  by  Gxu°s  the 
mapping  71r0(R2  \  iZ2)  ~ ^  C,  defined  by  the  formula 

G^(*)  =  Gu,5(*  +  xo).  (106) 

Obviously,  G^s  =  Gu,s,  and  Gx°s  will  also  be  referred  to  as  a  Modified  Gaussian  Beam. 

The  following  lemma  shows  that  for  large  u,  S,  the  values  Gu(x)  are  almost  zero  for  all 
x  €  r3.  It  is  an  immediate  consequence  of  Lemma  3.2  and  the  formulae  (103),  (104). 

Lemma  4.1  Suppose  that  u,S  are  two  positive  real  numbers,  such  that  6  <  uf  2.  Then  for  any 

x  e  r3, 

I  G.(x)  |<  (107) 
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The  following  lemma  shows  that  for  large  u,  6,  the  functions  GUts,  Gu  very  nearly  coincide 
for  values  of  argument  outside  IV 

Lemma  4.2  Suppose  that  u,6  are  two  real  positive  numbers,  such  that  6  +  1  <  u/2.  Then  for 
any  x  €  R2  \  Ra, 

|Gu,s(x)-Gu(i)|<e-fi.  (108) 


Proof. 

Since  Gu  is  a  solution  of  the  equation  (24),  analytic  outside  R\  and  satisfying  the  radiation 
condition  (26)  at  oo,  Theorem  2.1  yields 

G„(x)  =  -i  ■  f  (<?„(•  1)  •  |£((,x)  +  ^%)  ■  G(t,  x))dl,  (109) 

for  any  x  €  R2  \  R\.  Obviously  (see  Figure  2), 

Jr  (^«(<)  •  -Qfl  (<>  X)  +  \t) '  G(t,  x))dl  = 

fr(G.m '  If (<’ x) + ^w^(t> ' G(‘- *))rf' + 

l  (G.w  •  !§(<■*)  +  ■  G(«, *)M'-  (110) 

On  the  other  hand,  it  immediately  follows  from  the  combination  of  (31),  (32),  and  (107)  that 

1  /r,(G“(,)'lf(‘’I)+^W^(')'G(<’X))‘"  !<*-*,  (HI) 

and  we  obtain  (108)  by  combining  (111)  with  (110)  and  (105). 

□ 

The  following  lemma  shows  that  the  far-held  representation  of  the  modified  Gaussian  beam 
GU}$  is  almost  identical  to  that  of  the  Gaussian  beam  Gu.  Its  proof  is  similar  to  that  of 
Lemma  4.2,  and  is  omitted. 

Lemma  4.3  Suppose  that  u,6  are  two  real  positive  numbers ,  such  that  6  +  1  <  u/2 .  Then 

I  Fo(Gu,s)(9)  -  |<  e~s.  (112) 

for  all  6  €  [ — tt,  7t]. 


19 


4.2  Modulated  Gaussian  beams 
We  will  define  the  function  :  R2  —  C  by  the  formula 
K(z,  y)  = 

r 171  i  r2ir 

— — —  •  v/(- — )  •  /  Gu(x  -  m  •  cos(ri),  y  -  m  •  •  e'^dr],  (113) 

and  the  function  :  R2  — ►  C  by  the  formula 

KAx'  y )  = 

r  i  ft* 

— )•  /  Gu>s(x-m-cos(ri),y-Tn-sin(v))-et'm'r>dri  (114) 

(the  definitions  (113),  (114)  are  correct  and  stable  due  to  Lemma  2.6).  We  will  be  referring 
to  as  Modulated  Gaussian  Beam,  and  to  M™g  as  Modulated  Modified  Gaussian  Beam 
(MMGB).  Obviously,  M™  is  a  radiation  potential  outside  the  region  Sm(Bu )  (see  (62),  (63), 
and  Lemma  2.1,  and  M™g  is  a  radiation  potential  outside  Sm(R 2)  (see  (101)).  The  following 
lemma  supplies  the  far- field  representation  for  M™. 

Lemma  4.4  For  any  real  u  >  0  and  integer  m, 

F0(M™)(9)  =  e*'m -6  •  (115) 

for  all  8  £  [— 7r,  7r], 


Proof. 

Combining  (113)  with  (61)  and  (57),  we  have 
Fo(M?)(6)  = 


dmi™-) 


Jm(m ) 

Jm(m) 


J[—) .  ft  FoiG^'00^ r')'m'st'n(’?)))(^  .  j-^dr,  = 

V  8  *  7T  Jo 

p2'7T 

Fo(Gu)(8)  ■  /  e‘-^-(('"-‘:o*(')).rn'«n(f!)).(cos(9)1»in(e)))  .  _ 

U  70 

«m  n  *  \ 

Jm(rn ) 

Jro(G'u)(0)  •  ft  e'-^i^M-cosW+sini^  sinte))  .  ei  m  rt ^  - 

Jo 

-rA  ■  fi-sA')  •  -MCJW  •  r 

J fny'ffi)  o  *  7r  */ 0 


(116) 
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Now,  we  obtain  (115)  by  combining  (116)  with  (48)  and  (64). 


□ 


The  following  lemma  is  an  immediate  consequence  of  the  combination  of  Lemmas  4.4,  4.3. 
It  supplies  the  (approximate)  far-field  representation  of  M™s. 

Lemma  4.5  For  any  real  u  >  0,  integer  m,  and  0  <  6  <  u/2, 

|  Fo(M%s)(0)  -  •  e«-<«*W-i>  |<  e~s,  (117) 

for  all  8  6  [— 7r,7r]. 

For  an  arbitrary  integer  n  >  0,  real  u  >  0,  and  real  0  <  6  <  u/2,  we  will  define  the  function 
Ql :  R2  -  C  via  the  formula 

£  jicw,  (ns) 


and  the  function  Q*<s  :  R2  — »  C  via  the  formula 


(119) 


Obviously,  Q "  is  a  radiation  potential  outside  Sn(Bu )  and  Q"  $  is  a  radiation  potential  outside 

Sn(Ri)- 

The  following  two  theorems  describe  the  far-field  representations  of  Q™,  Q*s,  respectively. 
They  follow  immediately  from  Lemmas  4.4,  4.5,  respectively,  and  the  obvious  fact  that 


y-  i-m-e  _  s*n((n +  -|).fl) 

m=-n  sin(l) 

Theorem  4.6  For  an  arbitrary  integer  n  >  0  and  real  u  >  0, 


(120) 


W)  = 


-  +  2)  •  *)  .  ,u.(co*(*)-l) 


5zn(-) 


(121) 


for  all  6  €  [-7T,  jr]. 

Theorem  4.7  For  an  arbitrary  integer  n  >  0,  real  u  >  0,  and  real  0  <  S  <  u/2, 


FoiQlsW)  -  |<  e~* 


(122) 


for  all  6  €  [ — 7r,  7t]. 
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The  following  theorem  is  an  immediate  consequence  of  the  combination  of  Theorem  3.5 
with  (119),  (114),  (108).  It  shows  that  the  support  of  the  function  Q^  s  is  shaped  like  a  beam, 
whose  width  is  closely  related  to  that  of  the  Gaussian  beam  Gu. 


Theorem  4.8  Suppose  that  n  >  0  is  an  integer  number,  and  u,S,  are  two  real  numbers  such 
that  0  <  S  <  u/2.  Then 

\Qls{x,y)\<e-s  (123) 

in  either  of  the  following  two  (intersecting)  regions: 

1.  x  €  [2  •  ,/u, 2  •  u],  |  y  |>  2  •  v'TlO  •  <5)  •  \/u  +  n. 

2.  x  €  [2  •  u,  oo],  |  y  |>  ^  •  x  +  n. 

The  following  theorem  is  obtained  from  the  preceding  one  by  elementary  algebraic  manip¬ 
ulation. 


Theorem  4.9  Suppose  that  under  the  conditions  of  Theorem  4-8,  the  numbers  u,  n,  are  defined 
by  (11),  (18),  respectively,  and  that  in  addition,  u>2-6.  Then 

i<?;,«(*.sOi<  (124) 


in  either  of  the  following  two  (intersecting)  regions: 

L  lyl>8-r. 

2.  x  €  [53,00],  \y  |>  8-V2-6  -f. 


Theorem  4.7  provides  an  analytical  expression  for  the  far-field  representation  Fo(Q^g)  of 
the  potential  Q*s  .  The  following  theorem  provides  a  somewhat  less  detailed  description  of 
the  Fourier  series  of  Fq(Q JJ  s). 

Theorem  4.10  Suppose  that  under  the  conditions  of  Theorem  4-  7,  r  >  0  is  a  real  number, 
and  the  numbers  u,  8,  v,  n,  m  are  defined  by  the  formulae  (11),  (10),  (3),  (12),  (13).  Then 
for  all  integer  j  such  that  |  j  |<  2  •  v, 

I  (Fo(Qi,,))i  —  1  l<  e.  (125) 

Furthermore,  for  any  j  such  that  \  j  |>  m/2, 

I  (fo(?;,())>  l<  £•  (126) 


Proof. 

Denoting  by  s,t  the  functions  [ — 7r,  7t] 


sin((n  +  \)  -0) 
sin(^) 


C  defined  by  the  formulae 


(127) 
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t(0)  =  e«-(«*W-1>,  (128) 

we  observe  that  the  Fourier  series  of  s,  t  are  given  by  the  formulae  (120),  (49).  Due  to  (127),  (128), 
we  can  rewrite  (122)  in  the  form 

I  -«(*)•«(*)  |<£,  (129) 

from  which  it  immediately  follows  (due  to  the  convolution  theorem)  that 

ll(W;,<))i-(<*l)ill<e.  (130) 

Now,  the  conclusion  of  the  theorem  follows  immediately  from  the  combination  of  Lemma  2.8 
with  (130),  (120),  (49). 

□ 


4.3  Proof  of  Theorem  1.1 

In  this  subsection,  we  use  the  analytical  machinery  developed  in  the  preceding  ones  to  prove 
Theorem  1.1  of  Section  1.  Given  real  numbers  p,  r,  £  satisfying  the  conditions  (1),  (2),  we  will 
define  the  mapping  /?£e  :  [-7r,7r]  — >  C  by  the  formula 

=  QZ,siP  ■  cos(6),  p  ■  sin(0)),  (131) 

with  6  defined  by  (10),  v  defined  by  (3),  u  defined  by  (8),  and  n  defined  by  (13).  The  following 
four  lemmas  show  that  the  function  satisfies  the  conditions  of  Theorem  1.1,  and  is  very 
close  to  the  function  defined  in  Section  1. 

The  following  lemma  is  obtained  from  Theorem  4.8  by  elementary  algebraic  manipulation. 

Lemma  4.11  Suppose  that  under  the  conditions  of  Theorem  4-10,  the  function  j3£c  :  [—7 r,  tt]  — ► 
C  is  defined  by  (131).  Suppose  further  that  either 

p  €  (16  •  r,-^— :],  and  \  6  |>  — (132) 

4*0  p 

or 

r  r2  .  .  Ul  6t 

/)€  — r,oo],  and  |  0  |> - + - ,  (133) 

4*o  r  p 

or  both.  Then 


I  KM  l<  e. 


The  following  lemma  is  an  immediate  consequence  of  the  preceding  one. 


(134) 
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Lemma  4.12  Suppose  that  under  the  conditions  of  Theorem  4-10,  the  function  6£  £  is  defined 
by  (131).  Then 


I  |<  £• 

for  all  0  6  [ — tt,  tt],  such  that 


>  8^2(1) 

'  r  p 


(135) 


(136) 


The  following  lemma  is  an  immediate  consequence  of  the  combination  of  (59)  and  Theo¬ 
rem  4.10. 


"S 


Lemma  4.13  Suppose  that  under  the  conditions  of  Theorem  4-10,  the  function  is  defined 
by  (131).  Then  for  all  j  €  [-2  •  u,  2  •  i/], 

I  (fa);  -  Hj(p)  |<  £.  (137) 

Finally,  Lemma  4.14  below  is  easily  obtained  from  the  combination  of  Lemma  2.13  with  (131),  (122). 

Lemma  4.14  Suppose  that  under  the  conditions  of  of  Theorem  4-10,  th>e  functions  6£e, 
are  defined  by  (131),  (15),  respectively .  Then  for  all  6  £  [— 7r,7r], 

I  -  K,M  l< E-  <138) 

Corollary  4.15  Obviously,  Theorem  1.1  is  an  immediate  consequence  of  Lemmas  4-12,  4-13, 

4-U ■ 


5  Numerical  Considerations  and  Experiments 

5.1  A  numerically  more  attractive  procedure 

Theorem  1.1  of  Section  1  provides  a  construction  of  the  function  satisfying  the  conditions  1. 
-  3.  of  Section  1.  However,  the  function  6££  supplied  by  Theorem  1.1  is  in  no  sense  optimal,  and 
in  fact  has  been  chosen  so  as  to  simplify  the  proof  of  Theorem  1.1,  not  to  lead  to  numerically 
most  efficient  schemes.  The  following  construction  turns  out  to  provide  a  function  b^e  that  is 
considerably  more  attractive  numerically  than  that  provided  by  Theorem  1.1. 

Given  real  numbers  r,p,e  such  that  0  <  r  <  p/4,  and  t  >  0,  we  will  define  the  integer  number 
v  as  the  smallest  real  number  such  that 

Jj(r)  <  £.  (139) 

for  all  j  >  v.  We  will  define  the  integers  m,  n  by  the  formulae 

m  =  6  •  v,  (140) 
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and  a  real  number  u  by  the  formula 


e“*  •/*(«)  =  £.  (142) 

Finally,  we  will  define  the  function  /£>r  :  [— 7r,  7t]  — ►  C,  by  the  formula  (14),  and  the  function 
bpr  t  :  [ — 7r ,  7r]  — *■  C  by  the  formula  (15). 

Observation  5.1  Obviously ,  the  above  procedure  defines  a  function  very  similar  to  that  pro¬ 
vided  by  Theorem  LI,  as  is  obvious  from  Lemmas  2.5 ,  2.8 .  However ,  the  analogue  of  The¬ 
orem  1.1  for  the  construction  of  this  subsection  is  somewhat  subtle ;  the  proof  fragments  into 
a  large  number  of  cases ,  depending  on  the  relative  sizes  of  r,  p ,  and  e.  On  the  other  hand , 
once  the  function  bpe  is  obtained ,  it  is  quite  trivial  to  verify  numerically  that  it  satisfies  the 
conditions  L,  2.  of  Subsection  LI,  which  are  the  two  conditions  necessary  for  bpc  to  be  & 
translation  operator  (to  a  fixed  precision  e).  Furthermore,  our  numerical  experiments  indicate 
that  in  most  cases,  the  above  construction  works  better  than  that  provided  by  Theorem  LI,  in 
the  sense  that  the  coefficients  p,q  in  the  formulae  (6),  (7)  are  much  smaller. 

Observation  5.2  Clearly,  evaluating  the  function  f£iT  at  m  equispaced  points  on  the  interval 
[ — 7r,  7r]  is  an  order  m  procedure  (see  (14));  given  f€yT,  the  function  bpc  :  [ — 7r,  7t]  — ►  C  (defined 
in  (15))  can  be  evaluated  at  m  equispaced  nodes  via  the  FFT,  provided  that  m  is  a  product 
of  powers  of  small  prime  numbers.  Thus,  in  our  computations,  we  altered  slightly  the  defi¬ 
nition  (140).  Specifically,  we  defined  m  as  the  smallest  positive  integer  that  is  a  product  of 
powers  of  2,  3,  and  5,  and  is  greater  than  6*^.  With  this  modification,  the  function  bpe  can  be 
constructed  for  a  cost  of  the  order  m  •  log(m),  which  is  sufficiently  fast  for  most  applications. 

5.2  Numerical  illustrations 

In  this  subsection,  we  illustrate  the  behavior  of  the  functions  /?£e  with  several  figures,  as  follows. 

In  Figure  3,  we  illustrate  the  behavior  of  the  function  Gu  with  u  =  1000,  by  plotting  the  loci 
of  points  x  in  the  plane  where  |  Gu(x)  |=  e,  with  e  =  l.Oi?  —  3,  l.Oi?  -  7,  l.Oi?  -  12.  The 
beam-like  structure  of  Gu  is  quite  transparent  from  this  plot. 

In  Figure  4,  we  illustrate  the  behavior  of  the  modulated  Gaussian  beam  by  comparing 
it  to  the  behavior  of  a  Gaussian  beam  GUj$  with  the  same  parameters  u,6.  Specifically,  we 
plot  the  loci  of  points  where  |  Q^s  |=  exp(-6),  \  Guj  |=  exp(-S).  In  the  regime  we  chose 
for  this  illustration,  both  functions  behave  very  much  like  expanding  beams,  with  the  same 
angle  of  expansion.  The  distance  between  the  graphs  is  roughly  equal  to  n,  as  it  should  be 
(see  (113),  (114),  (118),  (119)). 

In  Figure  5,  we  illustrate  the  behavior  of  the  functions  /3£s  as  p  grows,  for  four  values  of  e 
(e  =  1.0 E  —  3,  1.015  —  6,  1.0.E  —  9,  1.0 E  —  12).  Specifically,  we  plot  the  size  (in  radians)  0 
of  the  region  around  0  where  |  |>  £,  viewed  as  a  function  of  p.  It  is  quite  easy  to  see  that 


the  region  shrinks  as  p  grows,  apparently  converging  to  some  constant,  depending  on  £.  This 
behavior  is  in  agreement  with  Lemma  4.12. 

In  Figure  6,  we  illustrate  the  behavior  of  the  functions  /?£e  as  a  tool  for  the  reduction  of  the 
computational  cost  of  the  Fast  Multipole  Method.  We  plot  the  ratio  of  the  number  of  nodes 
required  to  discretize  the  function  f3$e  at  the  Nyquist  frequency  to  the  number  of  nodes  required 
to  discretize  the  function  vn  of  [6],  to  the  same  accuracy  and  under  identical  conditions.  The 
ratio  is  plotted  for  £  =  1.0E  -  3,  1.0£  -  6,  1.0£  -  9,  1.0 E  -  12,  r  =  100,  and  p/r  €  [4,40].  It 
is  easy  to  see  that  under  these  conditions,  the  improvement  is  quite  dramatic  for  low-accuracy 
calculations.  When  the  desired  precision  is  high,  the  improvement  is  much  less  impressive. 

Figure  7  is  similar  to  Figure  6,  except  that  here,  r  =  500.  Obviously,  in  this  case,  the  reduc¬ 
tion  in  the  computational  cost  is  much  greater  than  for  r  =  100;  this  is  in  agreement  with 
Theorem  1.1. 

Figure  8  shows  the  plots  of  the  absolute  values  of  the  function  /J£e  :  [ — 7r,  7r]  — ►  C  with  r  =  100, 
e  =  1.0E  -  6,  and  p  =  400,  1000,  10000.  Here,  it  is  obvious  that  /?£e  is  structured  like  a  bell, 
with  the  width  of  the  bell  decreasing  as  p  increases.  By  the  time  p  ~  r2,  the  bell  shrinks  to  a 
point. 

Finally,  Figure  9  contains  a  plot  of  the  real  part  of  the  function  :  [— 7r,7r]  — ►  C.  The 
function  is  so  oscillatory  that  this  plot  is  not  very  informative.  However,  it  is  clear  from  Figures 
8,  9  that  while  the  absolute  value  of  /?£e  looks  very  much  like  a  bell,  the  function  itself  is 
quite  oscillatory,  except  near  6  =  0. 

6  Generalizations  and  Applications 

Obviously,  the  purpose  of  this  paper  is  purely  technical  -  to  construct  numerical  tools  to  be  used 
in  the  design  of  Fast  Multipole  Methods  for  the  Helmholtz  equation.  Furthermore,  in  a  vast 
majority  of  applications,  the  problems  are  three-dimensional,  so  that  the  principal  (though 
by  no  means  the  sole)  purpose  of  a  two-dimensional  scheme  is  to  serve  as  a  model  before 
three-dimensional  algorithms  are  attempted. 

The  construction  of  this  paper  is  trivially  generalized  to  arbitrary  real  Helmholtz  coefficients 
by  rescaling.  The  construction  extends  to  complex  Helmholtz  coefficients  easily,  as  long  as  the 
real  part  of  the  Helmholtz  coefficient  is  positive;  in  this  case,  the  proofs  have  to  be  modified 
slightly.  The  construction  becomes  numerically  unstable  for  Helmholtz  coefficients  with  large 
negative  imaginary  parts. 

Our  numerical  experiments  show  that  the  construction  of  the  preceding  section  can  be 
sharpened  somewhat,  especially  for  relatively  small  r  and  p.  In  other  words,  there  exist  versions 
of  the  function  that  have  the  same  frequency  content  as  those  constructed  in  the  preceding 
section,  and  that  are  small  on  a  greater  part  of  the  interval  [-tt,  jt].  However,  in  our  experiments 
we  used  optimization  techniques  to  construct  such  functions,  at  a  significant  cost  in  CPU  time 
(the  cost  of  our  procedure  is  of  the  order  0(m3),  with  a  fairly  large  constant).  At  the  present 
time,  the  possibility  of  more  efficient  schemes  of  this  type  is  under  investigation. 
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Extension  of  the  results  of  this  paper  to  three  dimensions  is  quite  straightforward,  and  a 
paper  reporting  it  is  in  preparation.  The  author  is  currently  in  the  process  of  incorporating 
the  construction  of  Section  5  of  this  paper  into  a  Fast  Multipole  scheme  for  the  solution  of 
two-dimensional  scattering  problems.  These  results,  and  their  extension  to  three  dimensions, 
will  be  reported  at  a  later  date. 

7  Acknowledgments 

The  author  would  like  to  thank  Professor  R.R.  Coifman  for  many  useful  discussions,  and  for 
his  continuing  interest  and  support.  He  is  indebted  to  Dr.  Matviyenko  and  Dr.  Wandzura  for 
their  help.  Special  thanks  are  due  Professor  W.C.  Chew,  Dr.  R.L.  Wagner,  and  Dr.  J.M.  Song 
for  providing  a  most  stimulating  competition. 


References 

[1]  M.  Abramovitz  and  I.  Stegun,  Handbook  of  Mathematical  Functions,  Applied  Math.  Series 
(National  Bureau  of  Standards),  Washington,  DC,  1964. 

[2]  R.  L.  Wagner,  W.  C.  Chew,  A  Ray-Propagation  Fast  Multipole  Algorithm,  Microwave  and 
Optical  Technology  Letters,  Vol.  7,  No.  10,  July  1994. 

[3]  J.M.  Song,  W.  C.  Chew,  Fast  Multipole  Solution  using  Parametric  Geometry,  Microwave 
and  Optical  Technology  Letters,  Vol.  7,  No.  16,  November  1994. 

[4]  R.  Coifman,  V.  Rokhlin,  S.  Wandzura,  Faster  Single-Stage  Multipole  Method  for  the  Wave 
Equation,  10-th  Annual  Review  of  Progress  in  Applied  Computational  Electromagnetics, 
Vol.  1,  pp.  19-24,  Monterey,  CA,  March,  1994,  Applied  Computational  Electromagnetics 
Society. 

[5]  N.S.  Koshliakov,  M.M.  Smirnov,  and  E.B.  Gliner,  Differential  Equations  of  Mathematical 
Physics,  North-Holland,  Amsterdam,  1964. 

[6]  V.  Rokhlin,  Rapid  Solution  of  Integral  Equations  of  Scattering  Theory  in  Two  Dimensions, 
Journal  of  Computational  Physics,  86(2)  :  414  (1990). 

[7]  V.  Rokhlin,  Diagonal  Forms  of  Translation  Operators  for  the  Helmholtz  Equation  in  Three 
Dimensions,  Applied  and  Computational  Harmonic  Analysis,  v.  1,  1993,  pp.  82-93. 

[8]  G.N.  Watson,  A  Treatise  on  the  Theory  of  Bessel  Functions,  Cambridge  University  Press, 
Cambridge,  1980. 


27 


0.30E+04 


O.OOE+OO 


A 


- 1 - - - - - 1 - . - . - 1 - 1 - . - 1 - 

O.OOE+OO  0.30E+04  0.60E+04  0.90E+04 


Figure  1:  Regions  A,  B,  with  u  =  10000,  6  =  40 
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Figure  2:  Rectangles  Ri,R.2,R3,R4\  drawn  not  to  scale 
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Figure  4:  Functions  Q”g,  G% 
S  =  13.82  ss  — i 


0.80E+00 


0.60E+00 


0.40E+00 


0.20E+00 


Figure  6:  Improvement  in  the  number  of  nodes  in  the  discretizations  of 
functions  /?£e,  as  a  function  of  p/r;  with  r  =  100,  and 
fi  =  1.0 E  -  3,  1.0 E  -  6,  1.0 E  -  9,  1.0E  -  12. 


Figure  7:  Improvement  in  the  number  of  nodes  in  the  discretizations  of 
functions  as  a  function  of  pjr\  with  r  =  500,  and 
€i  =  1.0 E  -  3,  l.OE  -  6,  1.0 E  -  9,  1.0 E  -  12 
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Figure  8:  Absolute  values  of  the  function  /?£e,  with  r  =  100, 
£  =  1.0E  -  6,  and  p  =  400,  1000,  ’lOOOO 


